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ABSTRACT 

Seismic  networks  often  tend  to  overestimate  the  magnitude  of 
earthquakes,  because  those  stations  within  the  network  that  do  not  detect  a 
particular  event  are  ignored  in  the  conventional  magnitude  averaging  procedure. 


By  assuming  a normal  distribution  of  world-wide  magnitudes  for 
any  given  event,  it  is  possible  to  establish  a simple  statistical  model  that  in- 
cludes the  additional  information  that  the  event  magnitude  at  non-detecting  sta- 
tions must  be  below  a certain  threshold  value. 

In  this  report,  maximum  likelihood  estimation  is  applied  to  de- 
termine event  magnitude  based  on  this  model.  The  advantages  and  limitations 
of  the  technique  are  discussed  using  both  simulated  and  real  data.  It  is  found 
that  the  maximum  likelihood  method,  when  applied  properly,  has  the  potential 
to  yield  a significant  improvement  in  network  magnitude  estimates  compared  to 
the  conventional  averaging  technique. 
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herein  which  has  been  supplied  by  other  organizations  or  contractors,  and  this 
document  is  subject  to  later  revision  as  may  be  necessary.  The  views  and  con- 
clusions presented  are  those  of  the  authors  and  should  not  be  interpreted  as 
necessarily  representing  the  official  policies,  either  expressed  or  implied,  of 
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SECTION  I 
INTRODUCTION 


When  estimating  the  magnitude  of  an  earthquake  recorded  by  a 
seismic  network,  the  common  approach  is  to  average  all  magnitudes  measured 
at  those  individual  stations  that  actually  detect  the  event.  This  procedure  often 
leads  to  overestimating  the  magnitudes  of  events  that  are  near  the  network  de- 
tection threshold  since  many  stations  will  not  detect  such  events,  and  there- 
fore be  ignored  in  the  averaging  procedure.  Clearly,  those  stations  will  us- 
ually be  the  ones  with  the  weakest  signal,  and  the  net  effect  is  to  introduce  a 
positive  bias  in  the  estimation  procedure. 

Herrin  and  Tucker  (1972)  computed  the  expected  error  intro- 
duced by  the  above  magnitude  estimation  method  for  the  case  of  a homogeneous 
or  near -homogeneous  network.  Their  basic  assumption  was  that  world-wide 

bodywave  magnitudes  of  a given  event  follow  a Gaussian  distribution  with  un- 

2 

known  mean  and  variance  (T  . They  called  ^ the  'true'  magnitude  of  the 
event,  and  computed  the  bias  relative  to  this  (unknown)  value  as  a function  of 
a2  and  the  network  characteristics. 

This  report  uses  the  same  basic  hypothesis  as  Herrin  and  Tuck- 
er regarding  the  normality  of  worldwide  event  magnitudes.  However,  the  esti- 
mation problem  is  approached  in  a different  way.  Specifically,  we  assume  that, 
for  a given  event,  the  magnitudes  at  all  detecting  stations  are  measured  as  well 
as  upper  limits  on  the  magnitudes  at  non-detecting  stations.  These  upper  limits 
are  obtained  by  measuring  the  largest  noise  peak  within  the  expected  signal  ar- 
rival window,  and  computing  the  corresponding  'noise  magnitude'.  Then,  infor- 
mally, we  find  the  'most  likely'  magnitude  based  upon  this  combined  informa- 
tion. More  precisely,  we  apply  the  statistical  method  of  maximum  likelihood 
estimation  to  the  above  problem. 
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Some  comments  are  appropriate  regarding  the  Gaussian  model 


for  world-wide  event  magnitude  distribution,  Freedman  (1967)  studied  the  am- 
plitude distribution  of  bodywave  recordings  published  by  the  United  States 
Coast  and  Geodetic  Survey  (USCGS),  and  found  a close  correspondence  with  the 
log  normal  distribution.  Clearly,  a lognormal  distribution  of  signal  amplitudes 
implies  a normal  distribution  of  magnitudes.  This  result  has  been  supported 
by  other  observational  studies.  The  normality  of  surface-wave  magnitudes 
has  not,  to  our  knowledge,  been  definitely  established,  but  still  seems  to  be  a 
reasonable  assumption. 

It  is  natural  to  ask  whether  this  statistical  picture  of  world-wide 
magnitude  distribution  is  consistent  with  recent  advances  in  earthquake  source 
theory,  or,  phrased  in  another  way,  if  most  of  the  observed  variance  could  be 
predicted  from  source  and  regional  information.  Clearly,  the  following  fac- 
tors influence  the  station  magnitudes  for  a given  event,  and,  if  known,  could 
be  used  to  compute  event  parameters  in  a deterministic  way  from  a limited 
set  of  observations: 

• Source  type  and  complexity  of  source  time  function 

• Near-source  earth  structure 

• Attenuation  of  signal  due  to  transmission  path  effects 

• Near-receiver  crustal  effects  such  as  scattering  of  P-waves. 

While  path  and  receiver  effects  may  be  accurately  predicted  if 
a large  number  of  events  from  a given  region  are  available  for  calibration  pur- 
poses, the  source  and  near-source  effects  on  the  amplitude  patterns  of  P-waves 
(in  particular)  are  very  difficult  to  specify.  For  example,  one  would  expect  an 
underground  explosion  to  produce  less  variation  in  world-wide  amplitudes  than 
an  earthquake;  still,  Lambert  et  al.  (1970)  found  the  Longshot  explosion  to  ex- 
hibit greater  amplitude  variation  than  the  average  earthquake. 
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It  is  also  difficult  to  infer  the  radiation  pattern  of  an  earthquake 
from  that  of  a previous  event  in  the  same  area.  For  example,  Ringdal  (1974) 
found  significant  scattering  in  the  magnitude  residuals  between  the  LASA  and 
NORSAR  arrays  for  a large  aftershock  sequence  from  the  Kurile  Islands  in 
June,  1973. 

Finally,  many  authors  now  believe  that  most  large  earthquakes 
are  composed  by  a number  of  events,  possibly  with  different  fault  planes  (Wyss 
and  Brune,  1967;  Blandford,  1974;  Burdick  and  Helmberger,  1974).  This  view- 
point implies  that  adequate  source  solutions  may  be  extremely  difficult  to  find. 
Clearly,  difficulties  also  arise  when  processing  smaller  earthquakes,  which 
may  be  simpler  in  nature,  but  lack  available  recordings  of  high  signal-to-noise 
ratio. 

In  summary,  the  statistical  approach  to  describing  the  world- 
wide event  magnitude  distribution  has  the  advantage  of  getting  around  the  con- 
siderable difficulties  inherent  in  computing  deterministic  event  radiation  pat- 
terns, while  still  retaining  a realistic  picture  of  actually  observed  data.  The 
true  magnitude  of  an  event  is  a parameter  of  the  assumed  statistical  distribu- 
tion, and  can  therefore  be  estimated  from  statistical  considerations  alone. 

Section  II  of  this  report  establishes  the  maximum  likelihood  ap- 
proach to  the  magnitude  estimation  problem.  In  Section  III,  simulated  event 
data  are  introduced  in  order  to  evaluate  the  estimator.  Section  IV  applies  the 
method  to  real  seismic  data  recorded  by  the  World-Wide  Standard  Seismograph 
Network  (WWSSN)  and  the  Very  Long  Period  Experiment  (VLPE)  network. 
Finally,  conclusions  and  recommendations  from  this  study  are  presented  in 
Section  V. 


SECTION  II 


DESCRIPTION  OF  THE  MAXIMUM  LIKELIHOOD  METHOD 


This  section  gives  a brief  description  of  the  statistical  model 
leading  to  the  development  of  the  maximum  likelihood  network  magnitude  esti- 
mator. The  likelihood  function  is  established  for  two  separate  cases,  corres- 
ponding to  whether  the  individual  station  detection  thresholds  are  known  by  ac- 
tual measurement  or  only  by  statistical  distribution.  The  asymptotic  proper- 
ties of  the  estimator  are  also  derived. 

A.  THE  GAUSSIAN  MODEL 

Our  basic  assumption  is  that,  for  a given  event,  world-wide 
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magnitudes  follow  a Gaussian  distribution  with  parameters  (j/.cr  ).  The  (un- 
known) mean  of  this  distribution,  /I  , is  called  the  true  magnitude  of  the  event, 
and  our  objective  is  to  estimate  this  parameter  based  upon  a set  of  observa- 
tions. 

Further,  we  assume  that  at  a given  station,  an  event  is  detect- 
ed if  the  station  magnitude  m exceeds  a certain  threshold  value  a . This 
threshold  value  may  be  treated  as  a random  variable,  if  desired,  but  we  will 
first  assume  that  a is  actually  measured  as  the  'noise  magnitude'  at  the  ex- 
pected time  of  signal  arrival. 

Thus,  we  assume  that  the  probability  of  detecting  a given  event 
may  be  written  as: 

P(  Detect/4,  a)  = P(m>a/ji,a)  = 0 (^~j  (II- 1) 

where  0 is  the  standard  cumulative  normal  distribution  function. 
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More  generally,  if  the  threshold  magnitude  m is  considered 

^ 2 

a normally  distributed  random  variable  with  mean  a and  variance  CT^,  , and 
if  m and  m,^  are  considered  independent,  we  have: 

P(  Detect/ , cr)  = P(m  > , a)  = 


= P(m-mT  > 0/n,  (j  ) 


(H-2) 


Since  m-m^,  is  a Gaussian  variable,  we  obtain 


P(Detect/^t , a ) = 0 


("if) 


(II- 3 ) 


where 


2 2 2 
cr  = ° + cr 

1 T 


(H-4) 


This  model  is  illustrated  in  Figure  II- 1, 

In  the  following,  we  will  concentrate  on  the  approach  with  a 
known  (non-random)  threshold  magnitude,  since  this  parameter  can  always  be 
measured.  Clearly,  it  is  more  satisfactory  to  know  the  precise  threshold 
magnitude  than  a statistical  distribution,  especially  when  taken  into  account 
that  the  statistical  distribution  may  not  always  be  valid;  e.g.  , in  cases  of  high 
coda  levels  following  a large  earthquake  or  in  the  event  of  instrument  malfunc- 
tioning. 


B.  THE  LIKELIHOOD  FUNCTION 

Assume  that  for  a given  event,  records  from  a network  of  n 
stations  are  examined.  Further,  assume  that  the  threshold  magnitudes  , 
i = 1,  2,...  n are  known,  and  that  for  those  stations  that  detect  the  event,  a 
magnitude  m.  is  computed.  Finally,  assume  that  all  station  observations 
may  be  considered  independent. 
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We  want  to  compute,  as  a function  of  the  unknown  event  para- 
meters (n,  a),  the  probability  that  a given  situation  occurs,  and  then  maxi- 
mize this  likelihood  function  with  respect  to  the  unknown  parameters.  A for- 
mal development  of  this  is  presented  in  Appendix  A;  the  likelihood  function  is: 


Mm  ...  m /jl.cr ) = 
i n 


where 


0(x)  = 


n -1-4— ) • n 'l 

All  ' ° I All  Non-\  ” ) 


<f>(x)  = 


Detections 


1 -x2/2 

y/z7  * e 


Detections 


■/ 


<t>( t)  dt  . 


(H-5) 


(n-6) 


(II—  7 ) 


It  is  of  interest  to  note  that  the  first  group  of  prr  iucts  in  (II- 5 ) 
represents  the  likelihood  function  in  the  special  case  that  all  non-detections 
are  ignored,  and  is  maximized  by  the  'conventional'  estimate  of  H ; i.  e.  , the 
average  of  the  magnitudes  at  the  detecting  stations. 

Since  the  second  factor  group  in  (II- 5)  is  a product  of  decreas- 
ing functions  of  /i  , it  follows  that  our  maximum  likelihood  estimate  will  al- 
ways be  less  than  or  equal  to  the  conventional  estimate,  with  equality  only  if 
all  stations  detect  the  event. 

Although  we  are  primarily  interested  in  the  parameter  , the 
likelihood  function  must  in  general  be  maximized  as  a function  of  two  paramet- 
ers (H,a)  . However,  it  may  in  some  cases  be  legitimate  to  restrict  a to  as- 
sume values  within  a certain  interval,  or  even  to  a predefined  value.  For  ex- 
ample, Veith  and  Clawson  (1972)  found  a value  of  cr=  0.4  to  be  representative 
of  the  WWSSN  short-period  network.  Bungum  and  Husebye  (1974)  showed  by 
comparing  NORSAR  and  NOAA  magnitudes  that  the  standard  deviation  of  0.  3 
appeared  to  be  independent  of  event  magnitude.  Clearly,  the  lower  standard 
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deviation  for  NORSAR  than  for  a WWSSN  station  is  reasonable  since  NORSAR 
magnitudes  are  averaged  over  a large  number  of  sensors  (more  than  100).  In 
conclusion,  it  appears  that  the  parameter  (X  may  be  predicted  reasonably  well, 
so  that  tne  variation  of  this  parameter  may  be  restricted  within  a fairly  nar- 
row range. 

For  an  inhomogeneous  network,  it  would  probably  not  be  right 
to  use  the  same  value  of  (X  for  each  station.  A more  general  approach  allow- 
ing individual  values  of  a as  well  as  individual  station  bias  values  is  present- 
ed in  Appendix  A. 


It  is  a straightforward  extension  of  the  preceding  approach  to 

find  the  likelihood  function  for  the  case  that  the  threshold  magnitudes  m 

2 Tl 

are  random  variables.  Thus,  if  m^.  is  N(a^,  <x  ) (i  = 1.2,...  n): 


JL^m^  . . 


where 


mjn.a)  - 


All  ° \ / All  Non-  \ ) 


( II- 8) 


All 
Detections 


Detections 


2 2,2 
(X.  = a + a . 
l Ti 


(II- 9) 


C.  APPROXIMATE  CONFIDENCE  LIMITS 

One  of  the  most  prominent  features  of  maximum  likelihood  esti- 
mators is  that  they  often  possess  very  desirable  asymptotic  properties.  Un  ler 
reasonably  general  conditions,  the  following  may  be  proved  (Cramer,  1945): 

• The  solution  of  the  likelihood  equation  converges  in  probability 
to  the  true  parameter  value  as  the  number  of  observations  in- 
creases 

• The  maximum  likelihood  estimator  is  asymptotically  efficient. 
(Informally,  an  efficient  estimator  is  one  that  has  a variance 
lower  than  any  other  unbiased  estimator) 
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• The  maximum  likelihood  estimator  is  asymptotically  normally 

distributed. 

In  order  to  verify  that  these  properties  apply  to  our  particular 
case,  we  find  it  convenient  to  regard  the  selection  of  network  stations  as  a 
random  experiment,  and  assume  that  the  probability  density  function  of  the 
threshold  magnitude  a is  of  some  fixed  form  s(a).  We  also  assume  that  in- 
dividual station  selections  are  independent. 

In  this  way,  we  can  view  the  estimation  procedure  as  consisting 
of  observing  the  outcomes  (m,  a)  of  n independent  experiments,  each  with  the 
likelihood  function: 


A(m,  a/ /i  , o ) = 


s(a)  • <J) 


(v) 

\ ° / 


for  m - a 


| s(a)  • 01 — • for  m£a  . 


(11-10) 


It  is  readily  verified  that  the  original  likelihood  function  (II- 5)  is  equivalent  to 
a product  of  n functions  of  the  form  (II-10),in  the  sense  that  the  factors  orig- 
inating from  s(a)  do  not  depend  upon  n and  a , and  therefore  will  not  in- 
fluence the  maximum  likelihood  estimates. 


As  n — oo , we  can  now  apply  the  two-dimensional  form  of  the 
limiting  theorem  in  Cramer  (1945),  in  order  to  show  that  the  asymptotic  pro- 
perties described  above  apply  to  our  case. 

The  important  implication  of  this  result  is  that  for  i given  event, 
adding  new  (independent)  stations  to  an  already  existing  network,  will  cause  the 
maximum  likelihood  estimate  to  converge  to  the  true  value,  even  if  the  new  sta- 
tions added  have  roughly  the  same  capabilities  as  the  original  network  stations 
(i.  e.  , s(a)  is  kept  constant).  An  example  of  this  principle  will  be  shown  in 
Section  III. 
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It  is  clearly  desirable  to  obtain  an  expression  for  the  confidence 
limits  of  the  maximum  likelihood  estimator.  In  view  of  the  preceding  consid- 
erations, it  is  reasonable  to  compute  the  variance  of  an  unbiased,  efficient 
estimator  of  /i  (usually  known  as  the  Cramer-Rao  bound)  and  use  this  value 
as  an  approximation.  Such  a computation  is  carried  out  in  Appendix  A,  under 
the  assumption  that  (T  is  known,  and  the  resulting  expression  is: 


Var  /i 


.)+(!-  0(z  ))  + 


[^<Zi,]‘ 


0(a.) 


(II- 11) 


where 


z . = 

1 


V/i 


i 1 , 2, . . . n 


(11-12) 


and  <p  and  0 are  defined  by  (II-6)  and  (II-7). 


Examples  of  using  this  expression  together  with  the  normal  dis- 
tribution to  approximate  the  distribution  of  the  maximum  likelihood  estimator 


will  be  studied  in  Section  III. 
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SECTION  III 

SIMULATED  PERFORMANCE  OF  THE  ESTIMATOR 


A.  INTRODUCTION 

In  general,  a closed-form  expression  of  a maximum  likelihood 
estimator  may  be  difficult  or  impossible  to  find.  Therefore,  the  exact  statis- 
tical distribution  of  the  estimator  usually  cannot  be  derived  (except  with  re- 
spect to  asymptot'c  properties).  In  many  cases,  the  most  practical  way  to 
determine  the  statistical  properties  of  the  estimator  is  to  simulate  its  perfor- 
mance in  selected  cases.  This  section  presents  simulation  results  for  the 
maximum  likelihood  network  magnitude  estimator  developed  in  Section  II. 

Several  different  simulation  experiments  were  conducted.  Ty- 
pically, the  procedure  for  each  of  these  was  as  follows: 

• Define  a hypothetical  seismic  network  with  known  threshold 

magnitudes  a.,  ...,  a for  each  of  the  n individual  stations 
1 n 

of  the  network 

• Select  an  event  magnitude  ^ and  standard  deviation  a , there- 
by assuming  that  the  distribution  of  actual  station  magnitudes  is 

known 

• Simulate  100  events  recorded  by  this  network.  For  each  event, 

n independent,  normally  distributed  random  numbers  x , x . 

12 

....  x^  are  generated  from  the  Gaussian  distribution  (/i,  cr^). 

These  numbers  are  assigned  as  station  magnitudes  for  the  par- 
ticular event 
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In  each  of  the  100  cases,  determine  detection/ no  detection  for 


each  of  the  n stations,  by  comparing  and  a (i=l,  2, . . . , n). 
Then  estimate  network  event  magnitude  in  the  conventional  way 
(by  averaging  over  all  detections)  and  by  maximizing  the  likeli- 
hood function  (II- 5 ) 

• Compare  the  resulting  100  conventional  and  maximum  likelihood 

estimates  to  the  expected  theoretical  distribution  of  event  mag- 
nitudes. 

In  most  of  our  simulation  experiments,  we  chose  a 10-station 
network,  with  threshold  magnitudes  in  even  increments  of  0.  1 from  4.  1 to  5.  0 
magnitude  units.  Thus  we  have  for  this  network,  which  we  call  Network  1: 

a.  = 4.  1 + (i-1)  * 0.  1 i = 1, 2, ...  10  (III- 1 ) 

This  is  thought  to  be  a realistic  approach  to  represent  the  detection  thresholds 
of  a reasonably  homogeneous  network  for  a given  seismic  event.  The  varia- 
tions in  threshold  values  would  in  practical  situations  be  attributed  to  differ- 
ences in  seismic  noise  levels,  epicentral  distances,  and  signal  radiation/pro- 
pagation effects. 

Another  network,  called  Network  2,  of  100  stations  was  also 
briefly  investigated.  The  threshold  levels  of  this  second  network  were  defined 
as: 


a = 4.1  + (j- 1 )/  1 0 * 0.  1 j = 1,2,...  100  (III-2) 

where  the  division  sign  denotes  integer  division.  Thus,  Network  2 had  10  sta- 
tions of  capability  4.  1,  10  stations  "f  4.2  etc.  up  to  10  stations  of  5.  0 thresh- 
old level. 
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Event  magnitudes  /i  were  selected  ranging  from  3.  5 to  5.  5. 

The  standard  deviation  o of  the  world-wide  magnitude  distribution  was  set  at 
0.4  magnitude  units  in  all  cases.  Finally,  we  applied  the  maximum  likelihood 
estimation  technique  both  with  o assumed  to  be  known  (Subsection  III-B)  and 
with  unknown  a (Subsection  III-C). 

B.  SIMULATION  WITH  KNOWN  STANDARD  DEVIATION 

The  first  group  of  simulation  experiments  was  carried  out  in 
order  to  determine  the  performance  of  the  maximum  likelihood  estimator  un- 
der ideal  circumstances,  i.  e.  , when  the  correct  value  of  the  standard  devia- 
tion (T  of  the  magnitude  distribution  was  known  a priori.  Thus,  the  likelihood 
function  (II-5)  in  these  cases  was  maximized  as  a function  of  one  variable  . 

The  results  are  presented  in  Figures  III- 1 through  III  - 5 , for 
five  different  event  magnitudes  H , ranging  from  5.  5 through  3.5.  Each  fig- 
ure covers  100  simulated  events  of  identical  (true)  magnitude,  recorded  by 
Network  1 (10  stations).  The  estimates  of  /i  resulting  from  applying  the  con- 
ventional method  are  shown  in  a histogram  at  the  upper  half  of  each  figure,  and 
the  corresponding  maximum  likelihood  estimates  are  shown  in  the  lower  half. 

The  two  smooth  curves  on  each  figure  are  Gaussian  distribu- 
tions centered  at  the  true  magnitude.  The  variance  of  the  upper  curve  is  the 
theoretical  variance  for  a magnitude  estimate  by  a network  of  n = 10  stations, 
provided  all  stations  detect.  The  variance  of  the  lower  curve  is  the  Cramer- 
Rao  bound  defined  by  equation  ( II - 1 1 ) , and  this  curve  therefore  gives  an  indi- 
cation of  how  well  the  maximum  likelihood  estimator  may  be  approximated  by 
its  asymptotic  limits. 

In  some  cases  an  event  was  not  detected  by  any  station  within 
the  network.  Such  events  are  marked  separately  on  the  figures,  and  of  course 
do  not  have  any  associated  magnitude. 
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CONVENTIONAL  ESTIMATES 

Average  Number  of  Detecting  Stations  = 9.  ?4  (97.  4%) 


EVENT  MAGNITUDE  DISTRIBUTION 
MAX  LIKE  ESTIMATES  (SI6MA=0,40I 

FIGURE  III-  1 

SIMULATED  PERFORMANCE  OF  NETWORK  1 FOR  100  EVENTS  OF 
MAGNITUDE  5.  5 AND  KNOWN  STANDARD  DEVIATION 
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Average  Number  of  Detecting  Stations  = 8.  14  (81.4%) 
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FIGURE  III  - 2 


SIMULATED  PERFORMANCE  OF  NETWORK  1 FOR  100  EVENTS  OF 
MAGNITUDE  5.  0 AND  KNOWN  STANDARD  DEVIATION 
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FIGURE  III  - 3 


SIMULATED  PERFORMANCE  OF  NETWORK  1 FOR  100  EVENTS  OF 
MAGNITUDE  4.  5 AND  KNOWN  STANDARD  DEVIATION 
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FIGURE  III-4 

SIMULATED  PERFORMANCE  OF  NETWORK  1 FOR  100  EVENTS  OF 
MAGNITUDE  4.  0 AND  KNOWN  STANDARD  DEVIATION 
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FIGURE  III- 5 

SIMULATED  PERFORMANCE  OF  NETWORK  1 FOR  100  EVENTS  OF 
MAGNITUDE  3.  5 AND  KNOWN  STANDARD  DEVIATION 
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The  major  results  from  these  simulation  experiments  are  as 
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• At  magnitude  5.  5 the  two  methods  are  essentially  equivalent 
and  unbiased  relative  to  the  true  magnitude 

• At  5.  0 and  lower  magnitudes,  the  conventional  estimates  ex- 
hibit a significant  positive  bias.  The  size  of  this  bias  shows  a 
gradual  increase  from  about  0.  1 magnitude  unit  at  5.  0 to  about 
0.  5 units  at  magnitude  4.  0 

• The  maximum  likelihood  estimates  are  clearly  superior  to  the 
conventional  estimates,  and  show  essentially  zero  bias  down 
to  magnitude  4.  0 


i 


• At  magnitude  3.  5,  both  methods  show  less  satisfactory  perfor- 

mance, mainly  due  to  the  low  detection  probability  at  this  mag- 
nitude 


• The  Cramer -Rao  bounds  appear  to  give  a good  approximation 

to  the  actual  variance  of  the  maximum  likelihood  estimator. 


In  order  to  get  an  impression  of  what  happens  when  the  number 
of  network  stations  increases,  while  the  average  individual  station  detection 
capability  remains  constant,  we  conducted  an  experiment  with  Network  2 (100 
stations)  as  shown  in  Figures  III- 6 and  III-7  for  M = 4.  0 and  3.  5,  respectively. 
Clearly,  the  number  of  non-detections  decreases  significantly.  More  interest- 
ing, the  following  principles  are  illustrated: 

• The  maximum  likelihood  estimates  converge  to  the  true  magni- 
tude as  the  number  of  stations  increases 

• The  conventional  estimates  converge  to  a n agnitude  value  that 
is  significantly  biased  relative  to  the  true  /alue  as  the  number 
of  stations  increases. 
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FIGURE  III  - 6 

SIMULATED  PERFORMANCE  OF  NETWORK  2 FOR  100  EVENTS  OF 
MAGNITUDE  4.  0 AND  KNOWN  STANDARD  DEVIATION 
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FIGURE  III-7 

SIMULATED  PERFORMANCE  OF  NETWORK  2 FOR  100  EVENTS  OF 
MAGNITUDE  3.  5 AND  KNOWN  STANDARD  DEVIATION 


Thus,  if  an  existing  network  is  augmented  with  new  stations  of 
about  the  same  detection  capability,  this  will  improve  the  maximum  likelihood 
magnitude  estimate,  but  not  reduce  the  bias  inherent  in  the  conventional  esti- 
mate. 

A final,  interesting  observation  is  that  it  is  possible  to  use  the 
maximum  likelihood  method  to  obtain  an  upper  bound  on  the  network  magnitude 
of  an  event  that  is  not  detected  by  any  individual  station.  This  may  be  achieved 
by  noting  that,  for  any  given  situation,  there  is  a certain  smallest  magnitude 
that  may  be  estimated  by  the  network.  This  magnitude  corresponds,  in  the 
case  of  either  Network  1 or  Network  2,  to  the  situation  where  one  station  de- 
tects the  event  with  magnitude  4.  1,  and  no  other  station  detects.  For  Network 
1,  the  maximum  likelihood  estimate  in  this  case  is  3.  8;  for  Network  2 it  is  3.4. 

Such  upper  bounds  on  non-detections  should  of  course  be  used 

with  caution,  like  all  statistical  estimates.  However,  properly  interpreted, 

they  may  be  of  value  when  considering  'negative  discriminants'  such  as  the 

absence  of  detectable  surface -wave s for  underground  explosions  in  M - m 

s b 

plots. 


C.  SIMULATION  WITH  UNKNOWN  STANDARD  DEVIATION 

The  likelihood  function  (II-5)  may  be  maximized  as  a joint  func- 
tion of  /!  and  a provided  that  the  network  consists  of  at  least  two  stations, 
and  that  at  least  one  station  detected  the  event.  (This  last  reuuirement,  of 
course,  always  applies.  ) This  subsection  presents  the  results  of  applying  two 
parameter  maximization  to  test  situations  analogous  to  thore  described  in  Sub- 
section III-B. 

It  is  intuitively  clear  that  any  estimate  of  O based  upon  only 
10  observations  would  usually  have  a fairly  large  error  margin.  This  is  par- 
ticularly pronounced  in  those  cases  when  only  a few  stations  detect.  In  order 
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to  minimize  the  effects  of  gross  error  in  the  estimate  of  a on  the  resulting 
estimate  of  the  magnitude  H , we  restricted  0 to  values  within  a predefined 
range  in  the  estimation  procedures.  This  range  was,  somewhat  arbitrarily, 
set  to 

0.  25  < a < 0.  60  . (Ill- 3) 

(We  recall  that  the  true  value  of  a is  0.40.  ) Restricting  a in  this  way 
should  not  make  a significant  difference  in  practical  applications  of  the  method, 
since  only  an  approximate  a priori  knowledge  of  signal  variance  is  required. 

Results  from  these  simulations  for  Network  1 are  presented  in 
Figures  III- 8 through  III- 1 1 for  event  magnitudes  U-  5.5,  5.  0,  4.  5,  and  4.  0, 
respectively.  The  following  major  points  may  be  made: 

• For  an  event  of  magnitude  5.5,  no  significant  difference  is  seen 
compared  to  the  case  of  known  o (Figure  III- 1 ).  This  is  con- 
sistent with  the  observation  that  the  maximum  likelihood  esti- 
mates of  H and  a are  independent  if  ali  stations  detect 

• For  H = 5.  0,  4.5,  and  4.0,  the  maximum  likelihood  estimates 
are  significantly  better  than  the  conventional  estimates.  Fur- 
thermore, it  appears  that  the  former  estimates  are  only  slight- 
ly inferior  to  those  made  with  a known  value  of  o’  , as  was 
shown  in  Figures  III -2  through  III-4 

• As  in  the  case  of  known  O , the  Cramer -Rao  bounds  appear  to 
give  an  adequate  picture  of  the  variance  of  the  maximum  likeli- 
hood estimator,  although  we  only  developed  those  bounds  for  the 
case  of  a known  o . 

As  a final  part  of  the  simulation  experiment,  we  investigated  the 
consequence  of  deliberately  using  a wrong  value  of  O when  maximizing  the 
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FIGURE  III- 8 

SIMULATED  PERFORMANCE  OF  NETWORK  1 FOR  100  EVENTS  OF 
MAGNITUDE  5.  5 AND  UNKNOWN  STANDARD  DEVIATION 
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FIGURE  III- 9 

SIMULATED  PERFORMANCE  OF  NETWORK  1 FOR  100  EVENTS  OF 
MAGNITUDE  5.  0 AND  UNKNOWN  STANDARD  DEVIATION 
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FIGURE  III- 10 

SIMULATED  PERFORMANCE  OF  NETWORK  1 FOR  100  EVENTS  C 
MAGNITUDE  4.  5 AND  UNKNOWN  STANDARD  DEVIATION 
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SIMULATED  PERFORMANCE  OF  NETWORK  1 FOR  100  EVENTS  OF 
MAGNITUDE  4.  0 AND  UNKNOWN  STANDARD  DEVIATION 


likelihood  function  (II-5).  The  effects  of  such  a mistake  will  be  most  pronoun- 
ced in  those  cases  when  few  stations  detect.  Figures  III- 1 2 and  III- 1 3 show 
the  resulting  estimates  for  M = 4.  0,  O = 0.  4 as  simulation  parameters,  with 
cr  set  to  0.25  and  0.60,  respectively  in  the  estimation  process.  A definite 
bias  is  seen  in  both  cases,  although  the  maximum  likelihood  estimates  are 
still  more  accurate  than  the  conventional  estimates. 

From  this  last  result  we  conclude  that,  unless  cr  is  known 
with  good  confidence,  the  best  way  in  practice  to  apply  the  maximum  likeli- 
hood estimation  method  is  to  use  a two-parameter  maximization  technique, 
and  allow  cr  to  vary  within  predefined,  reasonable  bounds. 
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SECTION  IV 
DATA  ANALYSIS 


In  this  section  the  maximum  likelihood  network  magnitude  esti- 
mation technique  is  applied  to  actual  seismic  data,  and  the  results  compared 

to  conventional  estimates.  Two  cases  are  discussed,  one  concerning  m esti- 

b 

mates  by  a subnet  of  the  World-Wide  Standard  Seismograph  Network  (WWSSN), 
the  other  relating  to  Rayleigh  wave  magnitudes  estimated  by  the  Very  Long 
Period  Experiment  (VLPE)  network. 

A.  APPLICATION  TO  A SUBSET  OF  WWSSN  STATIONS 

In  order  to  test  the  effectiveness  of  the  maximum  likelihood 
method  in  estimating  earthquake  bodywave  magnitudes,  the  natural  approach 
was  to  apply  the  method  to  data  from  the  World-Wide  Standard  Seismograph 
Network  (WWSSN).  Unlortunately,  a number  of  stations  in  this  network  do  not 
routinely  report  amplitude  and  period  with  each  detection.  Consequently  a very 
limited  number  of  magnitude  estimates  are  available,  and  stations  not  reporting 
magnitude  cannot  automatically  be  classified  as  not  detecting.  Furthermore, 
seismic  noise  levels  are  not  reported  by  these  stations,  and  must  be  estimated. 

The  following  approach  was  undertaken  in  order  to  minimize 
these  problems: 

• First  of  all  we  selected  as  a reference  event  set  an  earthquake 

aftershock  sequence  /from  the  Tadzhik-Sinkiang  border  region, 
August  11-31,  1974).  This  provided  a large  number  of  events 
(about  60)  from  the  same  seismic  region,  and  within  a relative- 
ly short  time  span,  during  which  the  WWSSN  could  be  assumed 
to  be  fairly  constant 
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Secondly,  a subset  of  13  WWSSN  stations  was  selected.  These 
stations,  which  are  listed  in  Table  IV-1,  were  the  ones  that 
appeared  to  report  the  largest  number  of  event  magnitudes  from 
the  aftershock  sequence.  The  large  arrays  NORSAR  and  L.ASA 
were  deliberately  not  included  in  this  subset 


For  each  of  the  above  13  stations,  we  assumed  a 'magnitude 
reporting  threshold1  to  be  the  average  value  of  the  three  lowest 
magnitudes  actually  reported.  This  provided  a set  of  numbers 
(Table  IV-1)  which  were  used  as  threshold  values  (a.,  i = 1,2, 

. . . , 13)  in  the  maximum  likelihood  estimation. 


We  thereby  obtained  a reasonably  homogeneous  network  with 
threshold  magnitudes  ranging  from  mb  = 4.  3 to  5.  3.  The  network  had  a bal- 
anced distribution  of  epicentral  distances,  and  all  stations  were  located  in  the 
tele  seismic  range  of  30-90  degrees.  The  azimuthal  distribution  was  not  quite 
so  good,  but  was  still  thought  to  be  adequate.  By  excluding  NORSAR  from  the 
network,  we  were  able  to  use  magnitudes  from  this  station  as  a reference  to 
check  the  validity  of  our  network  estimates. 

Table  IV-2  lists  the  events  of  the  reference  set.  All  events  re- 
ported in  the  PDE  bulletins  for  this  aftershock  sequence  within  the  time  inter - 
val  August  11-  August  31,  1974  are  included,  except  for  one  event  occurring 
while  NORSAR  was  out  of  operation  and  three  events  that  were  not  detected  by 
any  of  our  13  stations.  Individual  station  magnitudes  are  included  in  the  table, 
as  well  as  network  magnitudes  estimated  both  by  conventional  averaging  and  by 
the  maximum  likelihood  technique. 


For  the  maximum  likelihood  method,  equations  ( II- 8 ) and  (II- 9) 

were  applied,  with  a preset  to  0.  40  and  <r  = 0.  2 (i  = 1 , 2 13).  The 

choice  of  a is  consistent  with  the  results  by  Veith  and  Clawson  (1972),  while 
mainly  reflects  variations  in  seismic  noise  level,  and  can  therefore  rea- 
sonably be  set  to  a small  value  for  the  aftershock  sequence. 


, 

a 


IV-2 


ST  OF  SELECTED  WWSSN  STATIONS  AND  THEIR  LOCATION  RELATIVE 
TO  A SELECTED  EARTHQUAKE  AFTERSHOCK  SEQUENCE 


t- 

s£> 

r- 

oo 

00 

rO 

vO 

in 

m 

O' 

O' 

00 

V 

M4 

Tf4" 

in 

'i* 

Tt4’ 

f—i 

O 

CM 

o 

O 

n- 

o 

in 

U4 

oo 

o 

cm 

o 

o 

r- 

in 

i-H 

CM 

CM 

in 

in 

CM 

ro 

rO 

cO 

CM 

CM 

m 

m 

m 

oo 

r“< 

oo 

oo 

rO 

r- 

r- 

m 

m 

CM 

cm 

sO 

NO 

CM 

rO 

O' 

o 

o 

CM 

m 

vO 

oo 

'O 

m 

rO 

sO 

sO 

n- 

t*- 

t'- 

r- 

f- 

oo 

■ii4 

w 

w 

w 

W 

W 

£ 

£ 

£ 

w 

£ 

£ 

£ 

w 

in 

CM 

in 

m 

in 

in 

sO 

in 

o 

00 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

rO 

oo 

oo 

O' 

m 

r- 

00 

O' 

sO 

CM 

o 

CM 

m 

»-H 

m 

CM 

O' 

f-H 

o 

i-H 

p-H 

£ 

£ 

£ 

£ 

£ 

£ 

£ 

£ 

Cfl 

£ 

£ 

£ 

£ 

£ 

ID 

O' 

rO 

in 

CM 

rO 

o 

— 1 

nO 

m 

in 

n- 

oo 

o 

o 

vO 

sO 

• 

oo 

in 

o 

■■i4’ 

CM 

• 

• 

o 

sO 

ID 

in 

m 

sO 

vO 

CM 

s£> 

vO 

vO 

in 

sO 

^cQOScQS^UcQkCQ^h 


■— i('jrOT}<tn\Or-aOO'0-i(\jcO 


TABLE  IV -2 

ELECTED  WWSSN  STATION  MAGNITUDES  AND  NETWORK  MAGNITUDE 
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Plots  of  the  estimated  network  magnitudes  versus  the  NORSAR 
m^  values  are  shown  in  Figures  IV-1  and  IV-2  for  the  conventional  method 
and  the  maximum  likelihood  method,  respectively.  It  is  clear  from  these  plots 
that  the  latter  method  gives  by  far  the  best  correspondence  between  the  network 
and  NORSAR,  and  it  is  concluded  that  the  maximum  likelihood  method  produces 
more  reliable  magnitudes  than  the  conventional  method  for  this  data  set.  The 
bias  of  abouf  0.  1 units  in  Figure  IV-2  is  probably  caused  by  the  beamform- 
ing loss  at  NORSAR,  which  is  not  compensated  for  in  NORSAR  m.  estimates. 

b 

It  is  observed  from  Figure  IV-1  that  the  conventional  estimates 

seem  to  saturate  toward  a magnitude  of  about  5.  0 as  the  NORSAR  m,  value  de- 

fa 

creases.  This  is  precisely  what  could  be  expected  from  our  simulaHon  results 
in  Section  III.  The  maximum  likelihood  estimates  will  also  reach  a saturation 
level  for  any  given  network,  but  this  level  is  at  least  half  a magnitude  unit  low- 
er for  the  particular  situation  discussed  here. 

It  is  appropriate  to  add  thc_t  we  also  attempted  to  apply  the  maxi- 
mum likelihood  method  with  other  parameter  settings  such  as  = 0 for 

i = 1,2,..,,  n and  a two  parameter  maximization  of  the  likelihood  function. 

The  results  were  essentially  the  same  as  the  ones  presented,  and  therefore 
indicate  a robustness  of  the  maximum  likelihood  method  which  was  also  demon- 
strated in  the  simulation  experiments  in  Section  III. 

As  a final  remark,  we  note  that  some  of  the  largest  events  were 
not  measured  at  one  or  two  stations;  apparently  due  to  high  coda  level  remain- 
ing from  a previous  large  event.  Although  these  non-detections  did  not  sub- 
stantially impact  the  maximum  likelihood  magnitude  estimates,  it  is  neverthe- 
less clear  that  the  assumed  Gaussian  distribution  of  threshold  magnitudes  does 
not  adequately  represent  such  cases.  Therefore,  we  feel  that  it  is  very  much 
preferrable  to  actually  measure  the  threshold  magnitude  for  each  non-detection 
if  at  all  possible. 
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FIGURE  IV -1 

WWSSN  SUBNET  mb  VERSUS  NORSAR  m.  FOR  EVENTS  FROM 
AN  AFTERSHOCK  SEQUENCE  USING  CONVENTIONAL 
NETWORK  MAGNITUDE  ESTIMATES 
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FIGURE  IV -2 
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APPLICATION  TO  THE  VLPE  NETWORK 


Data  from  the  Very  Long  Period  Experiment  (VLPE)  network 
have  recently  been  extensively  analyzed  by  Lambert  et  al.  (1974).  Their  data 
base  consisted  of  more  than  1000  seismic  events  from  Eurasia  during  1972, 
1973,  and  1974.  The  11  individual  stations  of  the  VLPE  network  are  listed  in 
Table  IV-3.  None  of  these  stations  were  operational  throughout  the  full  time 
period  covered,  and,  most  often,  data  for  a given  event  were  available  only 
from  five  or  fewer  stations.  Hence  the  VLPE  network  was  not  ideal  for  our 
purpose,  but  was  still  found  useful  to  illustrate  some  aspects  of  the  maximum 
likelihood  estimation  technique. 

Our  procedure  in  processing  the  VLPE  data  was  as  follows: 

• Individual  VLPE  station  Rayleigh  wave  data  were  obtained  from 
all  events  in  Lambert's  data  base  that  were  not  classified  as 
presumed  explosions 

• For  each  event,  a subset  of  stations  was  selected  by  eliminating 
all  stations  that  either  were  non-operational,  had  poor  data 
quality  such  as  spikes,  or  were  influenced  by  the  coda  level  of 

a previous  large  event 

• The  resulting  sub-network  for  each  event  was  split  into  detect- 
ing and  non-detecting  stations.  For  each  detecting  station,  the 
Rayleigh  wave  magnitude  measured  by  Lambert  et  al.  was  as- 
signed. This  usually  was  the  20  second  period  M . If  no  20 
second  magnitude  had  been  measured,  Mg  at  30  or  40  seconds 
periods  were  used  instead 

• For  each  event  in  this  data  base,  Mg  values  were  estimated  by 
conventional  averaging  over  all  detecting  stations  and  also  by 
applying  maximum  likelihood  estimation 
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• In  the  maximum  likelihood  estimation,  the  approach  given  by 
(II-8)  and  (II- 9)  was  utilized.  A signal  standard  deviation  of 

<J  = 0.  40  was  used,  together  with  values  of  <7^.  = 0.  40  (i  = 1, 

• • • , 11).  These  last  values  are  consistent  with  observed 
variations  in  noise  level  at  the  VLPE  sites  (Prahl,  1974),  while 
the  value  of  (7  was  chosen  somewhat  arbitrarily 

• The  actual  detection  threshold  estimates  at  each  individual  sta- 
tion (a.)  at  60  degrees  distance  are  listed  in  Table  IV-3.  These 
values  were  obtained  from  Lambert  et  al.  , 1974,  using  the  50 
percent  detection  threshold  of  each  station  and  compensating 
for  the  average  epicentral  distance  to  the  event  set.  In  our 
maximum  likelihood  estimation  model,  the  actual  station  thresh- 
old a^  for  a given  event  was  found  as: 

ai  = ai  - log  60  + log  A (IV- 1) 

where  A is  the  epicentral  distance  from  the  station  to  the 
event  in  degrees. 

We  refer  to  Lambert  et  al.  (1974)  for  a comprehensive  list  of 
individual  station  data  and  event  parameters  for  the  VLPE  event  set. 

Figures  IV-3  and  IV-4  give  a comparison  between  the  results 
from  the  two  methods  of  estimation.  Each  of  these  figures  shows  the  seismic- 
ity line  and  the  incremental  detection  curve  for  the  event  data  set,  and  has  been 
derived  by  applying  the  method  of  Lacoss  and  Kelly  (1969).  The  seismicity  line 
is  defined  by  the  magnitude  frequency  relationship: 

log10Nc  = A “ B ’ M (IV-2) 

where  Nc  is  the  number  of  events  with  surface -wave  magnitude  exceeding  a 
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FIGURE  IV -3 

ESTIMATES  OF  SEISMICITY  AND  NETWORK  DETECTION  PROBABILITY 
(AT  LEAST  ONE  STATION  DETECTING)  FOR  THE  VLPE  NETWORK 
BASED  UPON  CONVENTIONAL  M ESTIMATES 
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FIGURE  IV -4 

ESTIMATES  OF  SEISMICITY  AND  NETWORK  DEFECTION  PROBABILITY 
{AT  LEAST  ONE  STATION  DETECTING)  FOR  THE  VLPE  NETWORK 
BASED  UPON  MAXIMUM  LIKELIHOOD  M ESTIMATES 


certain  value  M.  The  incremental  network  detection  curve  is  assumed  to  be 


a cumulative  Gaussian  distribution  function  of  the  form: 


(M-M  y 
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(IV-3) 


where  M is  the  50  percent  network  detection  threshold  and  s is  a standard 
o 


deviation. 


In  this  way,  Figures  IV-3  and  IV-4  are  interpreted  as  both  re- 
presenting the  estimated  seismicity  and  the  VLiPE  network  detection  capabil- 
ity based  on  the  same  set  of  data.  However,  in  Figure  IV-3  conventional  net- 
work magnitude  calculations  have  been  used,  while  Figure  IV-4  refers  to 
max. mum  likelihood  magnitudes.  Two  important  differences  are  observed  be- 


tween the  two  cases: 


The  maximum  likelihood  approach  yields  a small,  but  signif- 
icant, decrease  in  the  estimated  value  of  B (0.85  versus 


0.93) 


The  detection  thresholds  are  lower  (by  0.  2 to  0.  3 M units) 


when  maximum  likelihood  magnitudes  are  used. 


Both  of  these  differences  were  expected  from  theoretical  con- 
siderations. The  important  questions,  of  course,  is  which  of  the  two  methods 
provides  the  most  reliable  estimates. 


For  the  detection  thresholds,  we  may  compare  the  results  from 
Figures  IV-3  and  IV-4  to  those  obtained  independently  by  Lambert  et  al.  (1974). 
The  50  percent  estimate  is  the  most  reliable  (Ringdal,  1974)  and  will  be  used. 
Lambert  et  al.  found  a 50  percent  incremental  detection  threshold  for  VLPE 


of  4.2  m units  by  using  the  WWSSN  combined  with  LASA  and  NORSAR  as  a 
b 


I 
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reference.  This  corresponds  to  about  3.4  Mg  units,  using  their  empirically 
derived  conversion  formula: 

Mg  = mb  - 0.  80+  r (4.  2<mb<5.5)  (IV-4) 

where  r is  a random  Gaussian  variable  of  zero  mean.  Using  NORSAR  and 
ALPA  Mg  values  as  a reference,  and  correcting  for  transmission  path  effects, 
they  obtained  another,  independent  estimate  of  3.  2 for  the  50  percent  thresh- 
old. These  values  are  in  good  correspondence  with  the  value  of  3.  38  obtained 
by  the  maximum  likelihood  Mg  estimation,  while  the  value  of  3.64  yielded  by 
the  conventional  technique  appears  to  be  too  high. 

With  respect  to  the  difference  in  B -value,  it  is  not  easy  to  find 
a 'correct'  estimate  to  compare  against.  Most  recent  stuoies  have  concentrat- 
ed on  b-values  for  the  corresponding  magnitude -frequency  relationship  for 
bodywaves.  For  example,  Bungum  and  Husebye  (1974),  using  NORSAR  m 

b 

data,  estimated  regional  b-values.  For  the  three  regions  where  most  of  our 
events  were  concentrated  (central  Asia,  southern-eastern  Asia  and  Japan- 
Kuriles -Kamchatka)  their  estimates  were  0.  85,  0.  86,  and  0.  84,  respectively. 
This  is  in  surprisingly  good  (and  probably  somewhat  coincindental)  corres- 
pondence with  the  maximum  likelihood  result  of  0.  85  (Figure  IV-4).  It  is 
significantly  lower  than  the  slope  of  0.  93  in  Figure  IV-I,  estimated  by  the 
conventional  method. 

Some  precautions  are  clearly  necessary  here.  First  of  all, 
one  might  ask  whether  the  NORSAR  values  represent  the  'true'  slopes,  since 
they  are  based  upon  measurements  at  only  one  station.  However,  it  was 
shown  by  Ringdal  (1974,  Appendix  A),  that  under  the  assumption  of  a Gaussian 
magnitude  distribution  model,  with  constant  variance,  the  estimate  of  the 
slope  b by  one  station  will  be  unbiased  relative  to  the  true  value.  In  fact, 
for  precisely  the  reasons  discussed  in  this  report,  single  station  seismicity 
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estimates  may  be  more  reliable  than  estimates  based  upon  network  magni 
tudes. 
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Secondly,  it  is  not  clear  that  B = b ; i.  e.  , that  the  slopes  based 
on  Mg  and  are  identical.  However,  in  the  magnitude  range  where  the  re- 
lationship (IV-4)  is  valid,  this  can  be  shown  to  be  the  case. 

In  conclusion,  we  feel  that  it  has  been  demonstrated  that  the 
maximum  likelihood  network  magnitude  estimation  method  gives  more  reliable 
estimates  than  the  conventional  averaging  technique,  and,  if  used  properly, 
can  provide  significant  improvements  in  statistical  analysis  concerning  event 
detection,  discrimination,  and  seismicity. 
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SECTION  V 

CONCLUSIONS  AND  RECOMMENDATIONS 
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The  maximum  likelihood  technique  has  been  found  to  potentially 
yield  a significant  improvement  in  network  magnitude  estimates  of  small  and 
medium  size  events  compared  to  the  conventional  method  of  averaging  the  mag- 
nitudes of  all  detecting  stations.  Such  an  improvement  has  been  actually  ob- 
served both  for  simulated  networks,  for  a 13  station  subnet  of  the  WWSSN,  and 
for  the  11  station  VLPE  network.  From  our  investigations  it  appears  that  the 
maximum  likelihood  method  applied  to  a 1 0 station  network  will  provide  essen- 
tially unbiased  magnitude  estimates  at  about  1 magnitude  unit  lower  than  the 
conventional  technique  does. 


We  recommend  that  further  research  be  carried  out  in  order  to 
obtain  more  complete  data  on  the  maximum  likelihood  method  and  its  under- 
lying assumptions.  Specifically,  the  following  topics  are  suggested: 

• Further  verification  of  the  Gaussian  model  for  world-wide  seis- 
mic magnitude  distribution  for  a given  event;  especially  for  sur- 
face wave  magnitudes 

• More  precise  determination  of  the  standard  deviation  O in  the 
above  distribution,  and  its  possible  variation  with  source  func- 
tion, event  depth,  magnitude,  and  seismic  region 

• Comparison  of  <7  for  array  station  networks  and  <7  for  single 
station  networks 

• Actual  application  of  the  maximum  likelihood  technique  to  exist- 
ing and  planned  networks.  For  such  applications,  it  is  recom- 
mended that 
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All  threshold  values  of  non-detecting  stations  should  be 
actually  measured  for  each  event  detected  by  the  network 
A narrow,  but  realistic  range  of  a should  be  specified, 
and  the  likelihood  function  (II—  5 ) should  be  maximized  as 
a function  of  two  parameters  (n,o). 

It  is  strongly  recommended  that  in  any  future  operational  net- 
work all  threshold  magnitudes  for  non-detecting  stations  should  be  measured 
along  with  t^e  magnitudes  of  the  detecting  stations.  It  is  felt  that  realistic 
estimates  of  the  magnitudes  of  small  and  intermediate  events  will  ultimately 
have  to  take  all  of  this  information  into  account,  whether  or  not  the  techniques 
and  models  used  in  this  report  are  to  be  applied. 
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APPENDIX  A 

DERIVATION  OF  CRAMER-RAO  BOUNDS  FOR  THE 
MAXIMUM  LIKELIHOOD  ESTIMATOR 


This  appendix  presents  a detailed,  formal  derivation  of  the 
Cramer-Rao  bounds  on  the  variance  of  the  maximum  likelihood  estimator  de- 
veloped in  this  report.  Since  it  was  established  in  Section  II  that  the  estima- 
tor is  asymptotically  normally  distributed  with  minimum  variance,  these  bounds 
may  be  used  as  approximations  to  the  actual  variance  in  practical  cases,  assum- 
ing that  the  number  of  stations  in  the  network  is  sufficiently  large. 

We  use  in  this  appendix  a slightly  more  general  problem  for- 
mulation than  in  Section  II,  in  that  we  allow  for  individual  station  differences 
in  magnitude  bias  and  signal  variance.  The  statistical  test  situation  is  thus 
summarized  as  follows: 

• Seismic  records  are  examined  for  a total  of  n stations 

• For  station  i (i  = 1, 2, . . . , n),  it  is  assumed  that  the  signal 
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magnitude  M.  is  normally  distributed  M.~N(/i+  b.  ,cr.). 

® Each  station  has  a certain  detection  threshold  a measured  at 

i 

the  time  of  the  event;  i.  e.  , the  station  detects  the  event  if 
mi"ai  and  does  not  detect  if  m.<a.  . Here,  m.  denotes  the 
actual  value  of  the  random  variable  M.  , and  is  defined  re- 
gardless of  whether  or  not  it  actually  can  be  measured 

• M^  and  M^  are  independent  random  variables  for  all  pairs 

(i.  j)  I i 4 j. 
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In  the  following,  we  assume  that  the  values  (b.  i = 1,2, 


, , . , n are  known  a priori,  and  that 


Threshold  magnitudes  a^ a^  are  actually  measured  for 

all  stations  at  the  time  of  expected  signal  arrival 

Actual  event  magnitudes  {m^;  i€*Z)}  are  measured  for  all  the 
stations  that  detect.  Here,  &D  denotes  the  set  of  indices  cor- 
responding to  the  detecting  stations. 


We  thus  assume: 


m.  't  a. 


for  all  i € 


m.  < a. 
1 1 


for  all 
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The  likelihood  function  L,  i.  e.  , the  probability  of  a certain  occurrence  of  de- 
tections, no  detections,  and  measured  values  m.  can  now  be  derived: 


L(m, , , , . , m la)  • dm,  • • • dm,  • • = 
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Prob(  O M.  € < m.  , m.  + dm.  >0  O M.  < a.  ) = 
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assuming  the  boundary  conditions  (A-l)  are  satisfied. 


Hence: 


Lfnij rrW a ) ~ 


H 0.(m.)  • JT  0 .(a.)  if  (A-l)  is  satisfied 
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elsewhere 
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where 


(m.-in+bj) 


^(m.)  = 


2tt«  cr. 
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( A-5) 


Some  remarks  about  the  likelihood  function  are  appropriate  at 

this  point: 

• Although  all  magnitudes  m, , . . . m and  threshold  values 

In 

al  * * * &n  are  *,nc^u{^ec^  in  the  boundary  conditions  (A-l)  of  the 
likelihood  function  L,  only  the  subsets  { m.  , it  2?  \ and 
{a  . actually  enter  the  expression  of  the  value  of  L 

( A-3).  Therefore,  JL,  as  a function  of  )i  , may  be  maximized 
as  long  as  the  magnitudes  of  the  detecting  stations  and  the 
threshold  values  of  the  non-detecting  stations  are  known 

• It  is  clear  that  (A-3)  also  applies  in  the  two  special  cases  when 
either  all  stations  detect  or  no  stations  detect.  Note  that  in  the 
latter  case,  the  likelihood  function  will  be  monotonically  de- 
creasing, and  hence  cannot  be  maximized. 

We  now  proceed  to  compute  the  Cramer-Rao  bounds  on  the 
variance  of  an  unbiased  estimator  /u  of  the  parameter  n . From  Cramer 
(1945),  the  minimum  variance  of  j i is  given  by 


( A-6) 


where  E denotes  the  statistical  expectation,  with  the  sample  space  consisting 

of  all  possible  combinations  of  station  magnitudes  m, , . . . m . Of  course 

In 
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with  all  threshold  magnitudes 

a particular  set  of  m 

1 n 

stations  count  as  detections. 


ai’**‘*  an  known  a priori,  a specification  of 
will  automatically  imply  knowledge  of  which 


In  order  to  evaluate  (A-6), 
set  of  functions  L.  given  by 


we  find  it  convenient  to  introduce  a 


L.(m./W  ) = / 

( *i(ai> 


if  m.  > a 
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if  m.  < a. 
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It  is  clear  from  the  definition  (A-3)  and  the  boundary  conditions  (A-l)  that  we 
may  write 


•Mm  » 


n 

mn^>  = IT  . 


(A-8) 


Furthermore,  since  L.  and  L.  for  i 4 j relate  to  two  independent  observa- 
tions, (A-8)  rnay  be  interpreted  as  a product  of  independent  functions  of  ran- 
dom variables  (strictly  speaking;  with  m.  replaced  with  M.).  The  log  likeli- 
hood function  becomes: 


Hence 


n 


Now,  it  is  true  in  general  that  for  a set  | 
variables,  we  have 
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Therefore: 
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a log  l 
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We  now  evaluate  the  relevant  terms,  and  first  observe  that 
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for  m.  < a. 
l i 


i l 


Thus: 


/ 


a log  l. 

• </>  ,(m.)  dm.  - 

OH  ill 


00 


■/ 


ai 


(f>  (m,)dm,  + 0,(a.) 

aM  i i i i i 


-1  a 


OU 


<P 


.(a.)  . / rf>.(m.)  dm. 

* i J 1 1 1 


-2-  ( 1 -$.(a.))  + -4r  <Ma.)  = o 

aii  i i a^  i i 


( A- 14) 


where  we  have  used  the  definition  (A-5)  of  <#>.  and  reversed  the  order  of  de- 
rivation and  integration.  It  is  interesting  to  note  that  (A-14)  is  valid  regard- 
less of  the  actual  type  of  probability  density  function  </>.  . 


A-5 


ing: 


The  remaining  terms  of  equation  (A-12)  are  found  by  evaluat- 


C^y  ■ fpgf 


• (f>.( m.)  dm. 
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By  inserting  the  actual  normal  probability  density  function  de- 
fined by  ( A-4)  in  (A-1 3)  we  obtain: 


d log  L. 
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l i 


(A- 1 6) 


for  in.  < a. 
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When  these  expressions  are  substituted  in  (A-15)  and  the  rom 
putations  carried  out,  we  arrive  at  the  following  result: 
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We  can  relate  this  expression  to  the  standard  Gaussian  distri- 
bution by  substituting 


Z. 
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Consequently: 
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where  <t>  and  <t>  represent  the  standard  (0,  1)  Gaussian  distribution  and  dei 
sity  functions,  respectively.  With  this  notation,  (A- 17)  becomes: 
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Further,  by  recalling  (A-6)  and  (A-12),  and  observing  that,  by  (A-14),  the 
cross  terms  in  (A-12)  vanish,  we  have  the  result: 


Var(M)  = 1 


/i 


d log  L. 


(a-2; 


where  the  individual  terms  may  be  computed  using  (A-21). 

The  following  comments  to  this  result  are  appropriate:  First 
let  a.  -»  -oo  , i = 1,2, ....  n;  this  is  equivalent  to  saying  that  all  stations  de 
tect  with  probability  1.  In  this  particular  case,  we  obtain  from  (A-22) 


Var  n = 


• i 0 • 

1=1  l 


( A-23 


Note  that  the  likelihood  function  (A-3)  is  simplified  in  this  case  by  allowing  m 

non-detections,  and  the  maximum  likelihood  estimator  can  be  computed  ana- 
lytically: 


I his  is  indeed  the  minimum  variance  unbiased  estimator  in  this  special  case, 
and  its  variance  is  given  by  (A-2  3).  A more  familiar  expression  is  found  by 


setting  all  b.  = 0,  i = 1 , 2 • • • n and  assuming  all  cr.  = <7  i = 1 , 2 • • • n.  We 


then  obtain  the  standard  expression 
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i=  1 


Var  M = — 
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for  estimating  the  mean  jLZ  of  a normal  distribution. 


However,  it  is  important  to  notice  that  in  a practical  situation 
(when  a.  ^ “ 80  I i = 1,2,...  n),  the  detection  threshold  of  the  stations  will  af- 
fect the  variance  of  the  estimate  of  \i  , even  for  a station  that  detects  a given 
event.  This  is  not  unreasonable,  since  the  variance  reflects  what  would  hap- 
pen if  (hypothetically)  a given  random  experiment  was  repeated  a large  number 
of  times.  Depending  upon  its  detection  threshold,  a station  would  then  be  ex- 
pected to  have  a certain  percentage  of  non-detections,  and  this  is  reflected  in 
the  expression  for  the  Cramer-Rao  bound. 


A quantitative  estimate  of  the  contribution  of  a given  station  to 
the  reduction  in  variance  of  the  estimate  of  ^ can  be  obtained  by  inspecting 
Table  A- 1 . It  is  seen  that  the  weighting  factor  in  (A-21)  is  close  to  1 unless 
the  station's  noise  level  is  significantly  higher  than  the  actual  estimated  event 
magnitude. 


TABLE  A- 1 


TABLE  SHOWING  VALUES  OF  THE  WEIGHTING  FACTORS  ENTERING 
INTO  THE  CRAMER-RAO  BOUNDS  ON  THE  VARIANCE  OF  THE 
MAXIMUM  LIKELIHOOD  ESTIMATOR 
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a - (b  +fl) 
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$ and  $ are  the  standard  Gaussian  distribution 
and  density  functions,  respectively. 


